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A computation of the dynamical structure factor of topologically dis- 
ordered systems, where the disorder can be described in terms of euclidean 
random matrices, is presented. Among others, structural glasses and super- 
cooled liquids belong to that class of systems. The computation describes 
their relevant spectral features in the region of the high frequency sound. 
The analytical results are tested with numerical simulations and are found to 
be in very good agreement with them. Our results may explain the findings 
of inelastic X-ray scattering experiments in various glassy systems. 

I. INTRODUCTION 

Inelastic X-ray scattering (IXS) experimentsEHi and inelastic neutron scatteringJl^O on struc- 
tural glasses and supercooled liquids provided useful information on the dynamics of their 
amorphous structure, at frequencies larger than 0.1 THz. Those experiments show a regime, 
when the wavelength of the plane wave is comparable with the inter-particle distance, where 
the vibrational spectrum can be understood in terms of propagation of quasi-elastic sound 
waves, the so-called high frequency sound. This high-frequency sound has also been observed 
in molecular dynamical simulations of strongEUll and fragile^llli liquids, and it displays sev- 
eral rather universal featuresiHiSBiiHl! (see however Ref.0 for a dissenting view regarding 
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the silica properties). In particular, a peak is observed in the dynamical structure factor at a 
frequency which depends linerly on the exchanged momentum p, in the region O.lpo — 1-QPoj 
Po being the position of the first maximum in the static structure factor. When extrapolated 
to zero momentum, this linear dispersion relation yields the macroscopic speed of sound. 
The width of the spectral line, T is well fitted by 

T(p) = Ap x , i«2, (1) 

with A displaying a very mild (if any) temperature dependence. Although the same scaling 
of the spectral line is found in hydrodynamics^, in this computation the proportionality 
constant is basically the viscosity, which shows a very strong temperature dependence which 
is not observed. Moreover, the same scaling of T has been found in harmonic Lennard- Jones 
glasses^, and one can safely conclude that the p 2 broadening of the high-frequency sound 
has a physical origin different from hydrodynamics. Other interesting problems related with 
the high-frequency vibrational excitations of these topologically disordered systems^ regard 
the origin of the Boson peak or the importance of localization properties to understand the 
dynamics of supercooled liquidsil. 

The variety of materials in which the p 2 broadening appears suggests a straightforward 
physical motivation. However, the simplest conceivable approximation, a wave propagat- 
ing on an elastic medium in the presence of random scatterers, yields Rayleigh dispersion: 
r oc p A . This result is very robust: as soon as one assumes the presence of an underlying 
medium on which the sound waves would propagate undisturbed, as in the disordered-solid 
model§!@, the p A scaling appears even if one studies the interaction with the scatterers 
non-perturbativelyill. In this paper we want to show that when the distinction between the 
propagating medium and the scatterers is meaningless (as it happens for topologically disor- 
dered systems), the p 2 scaling is recovered. Recently, it has been shown0 that the numerical 
solution of the mode coupling equations, modified for the study of the glassy phase, describes 
in a qualitatively correct way the range of frequencies explored by IXS scattering, including 
the high frequency sound and the boson peak. In that theoretical approach though, one 



cannot find a clear scaling law, T being proportional to p 2 only at very low momentum and 
following a more complicated law for higher values. 

We want to investigate the problem from the point of view of statistical mechanics of 
random matricesil, by assuming that vibrations are the only motions allowed in the system. 
The formalism we shall introduce, however, is not limited to the investigation of the high 
frequency sound and it could be straightforwardly applied in different physical contexts. 

Let us look more carefully at the relation between the high frequency behaviour and 
vibrational dynamics in glasses and at the relation between vibrations and random matrices. 
The dynamical structure factor for a system of N identical particles is defined asB 



where (. . .) denotes the average over the particles positions r 3 - with the canonical ensemble. 

A well known approach to the calculation of the S(p,u) in the high-frequency (THz) 
region, consists in taking into account only the vibrational modes of the system. The basic 
assumption of this calculational strategy is that there is a clear separation among the time 
scales of the 'fast' degrees of freedom, which are supposedly modeled as vibrations, and the 
'slow' ones, which give rise to the diffusion of the particles. One expects this assumption to 
be reasonably accurate in the glass and in the super-cooled liquid phase. 

When focusing only on the short time region, the slow motions can be thought as 'frozen', 
hence one studies the displacements u around 'quenched' positions x, by writing the position 
of the z-th particle as rj(t) = X{ +v,i(t), and linearizing the equations of motion. Then one is 
naturally lead to consider the spectrum of the Hessian matrix of the potential, evaluated on 
the positions x of the particles. Calling uj 2 , the eigenvalues of the Hessian matrix and e n (i) 
the corresponding eigenvectors, the one excitation approximation to the S(p,u) at non zero 
frequency is given in the classical limit by: 




(2) 




n=l 



(3) 
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However, in the supercooled liquid phase, one cannot always assume that in a typical 
configuration at equilibrium all the normal modes have positive eigenvalues, negative eigen- 
values representing the situation where some particles are moving away from the position 
x. This observation has motivated the so-called instantaneous normal modes (INM) ap- 
proachii^ii in the study of supercooled liquids. In that approach the normal modes in a 
typical configuration are supposed to describe the short time dynamics of the particles in 
the liquid phase, while it has been suggested^ that diffusion properties could be studied by 
considering the localization properties of the normal modes of negative eigenvalues. 

On the other hand, it has been argued that in a very broad class of glassy systems 
below a critical temper at ureil, called dynamical temperarture Tjj (for spin glasses) or mode 
coupling temperature T MC t (for structural glasses), the particles spend the most of their 
time around the minima of the potential energy surface. Recent numerical simulations seem 
to confirm that picture0H. 

As a consequence, in the glassy phase the high frequency window (> .1 Thz) of the 
dynamical structure factor is supposed to be well described by the propagation of quasi- 
harmonic excitations around a quenched structure^. By assuming that this structure cor- 
responds to one of the many minima of the potential energy, one can introduce a harmonic 
approximation where only the vibrations around these minima are considered, and all the 
dynamical information is encoded in the spectral properties of the Hessian matrix evalu- 
tated on the rest positions, where it has not negative eigenvalues. It has been shown using 
molecular dynamics simulation that below the experimental glass transition temperature the 
thermodynamical properties of typical strong glasses@ are in a good agreement with such 
an assumption. Let us note that, within this assumption, the dynamical structure factor 
(||) depend on the minimum around which the particles oscillate, hence the knowledge of 
an infinite number of disordered positions x seems to be needed. Here we shall make the 
reasonable assumption that the dynamical structure factor is self-averaging, a type of as- 
sumption which generally turns out to be correct for macroscopic observables in disordered 
systems. 
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Hence we have turned our problem into the study of spectral properties of a kind of 
random matrix, once the probability distribution of the elements of the matrix is given. 
More precisely, the problem of the high-frequency dynamics of the system can be reduced, 
in its simplest version, to the consideration of a very particular type of random matrices, 
the euclidean matrices^, in which the entries are deterministic functions (the derivatives 
of the potential) of the random positions of the particles. If the system has a conservation 
law such as momentum conservation in our case, due to translational invariance, all the 
realizations of the random matrix have a common normal mode with zero eigenvalue: the 
uniform translation of the system. An euclidean matrix is determined by the particular 
deterministic function that we are considering, and by the probabilistic distribution function 
of the particles. However, for the sake of simplicity we shall concentrate here on the simplest 
kind of euclidean matrices!^. We consider N particles placed at positions Xi,i = 1,2, ...,N 
inside a box of volume V, where periodic boundary conditions are applied. Then, the scalar 
euclidean matrices are 

N 

M ij = 8 ij ^2f{x i -x k )-f{x i -x j ) , i,j = 1,2... ,N . (5) 
k=i 

In the above expression, f(x) is an arbitrary but deterministic scalar function depending on 
the distance between pairs of particles, and the positions {x} of the particles are drawn by 
a probabilistic law P[x}. Notice that the matrix (||) preserves translation invariance, since 
the uniform vector eo(i) = const is an eigenvector of M with zero eigenvalue. Let us note 
furthermore that there are not internal indices, since for the sake of simplicity the particle 
displacements are restricted to be all collinear. Of course, in the study of the vibrations 
of glasses, one should take into account the vectorial nature of those vibrations and split 
the Hessian matrix in its transversal and longitudinal (with respect to the directions of the 
external momentum p) parts. A careful analysis of that feature will be presented elsewhere^. 

Simple as these matrices are, many physical problems ranging from vibrational spectra in 
amorphous systems and electron hopping in amorphous semi conductor s@, to instantaneous 
normal modes (INM) in liquidsiHil and combinatorial optimization^, can be described with 



them using an appropriate choice of P. In this paper we will study the simplest case, in 
which P is a uniform distribution function. 

The dynamical structure factor for a scalar euclidean-matrix is given by 



S E {p,E) = Y,QnW{E-E n ), (6) 



Quip) = j-f 



N 



i=i 



(7) 



S(p,u) = 2ujS E (p,u 2 ), (8) 

where the overline stands for the average over the particles position and we have given the 
definition either in the eigenvalue space (Se(p,E)) and in the frequency space (S(p, u)). 

Finally, let us mention that not all the topologically disordered systems are represented 
by an euclidean random matrix. For instance, diffusion on random graphs and on small- 
world network@ where the connecting sites are chosen with a probabilistic law, are described 
by the Laplacian operator Wij . The matrix element Wij , % ^ j is one if i and j are connected 
and otherwise, while the diagonal part is tailored to have DjWy = 0, as required by the 
conservation of the number of diffused particles. 

The layout of the rest of this paper is organized as follows. In section we shall apply an 
expansion in the inverse of the particle density formulation to the problem of calculating the 
Se(p, E). The same calculation can be performed using the field theory of Ref.0, adapted to 
this problem^, with identical results^. The analytical computations will be confronted with 
a numerical simulation in section [ITJ. In section [TV] we shall briefly consider the problem 
of including particle correlations on our calculation, which shall produce some qualitatively 
new results, regarding the exponential tail of the density of states. In section we shall 
present our conclusions. 

II. THE 1/p EXPANSION 

The basic object that we will calculate is the resolvent 
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G(p,z) = — J^eiP^-Xj) 



1 



Yz-M 



(9) 



from which the dynamical structure factor is recovered using the distribution identity 
l/(x + z0+) = P{l/x) - m8(x): 



Se(v, E) = lim Im G(p, E + irj) . 

71 r)^0+ 



The resolvent will also give us access to the density of states of the system: 

1 



(10) 



iV 



9e(E) 



lim Im 2J 



i=l L 



E + ir) — M 



— lim Im lim G(p, E + in) . 

7T v^0+ p^oo 



(11) 

(12) 



We shall set up a perturbative expansion in the inverse density of particles, 1/p. The 
zero th order of this expansion (to be calculated in subsection |11 A| ) corresponds to the limit 
of infinite density, in which the system is equivalent to an elastic medium. In this limit the 
resolvent turns out to be extremely simple: 
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G(p,z 



z - e(p) ' 
<P) = p[/(0) - /(p)] • 



(13) 
(14) 



In the above expression f(p) is the Fourier transform of the function / defined in Eq.(^), 
that due to its spherical symmetry, behaves like e(p) = 0(p 2 ) for small p. We see that the 



dynamical structure function has two delta functions at frequencies linear in p: uo = ±y e(p). 
It is then clear that Eq.flTB]) represents the undamped propagation of sound in our elastic 
medium, with a dispersion relation controlled by the function e(p) . The order 1/p corrections 
to Eq.flTB]) will be calculated in subsection |11B| . They take the form of a complex self-energy, 
£(p,s), yielding 

1 



G(p, z) 
S E (p,E) 



z - e(p) - E(p, z) ' 
1 ImE(p,E) 



tt(E- e(p) - i?eS(p, E)f + (JmS(p, E)f ' 



(15) 
(16) 
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The dynamical structure factor, is no longer a delta function, but close to its maxima it has 
a Lorentzian shape. From Eq.(|i6|) we see that the real part of the self-energy renormalizes 
the dispersion relation, and thus it should behaves proportionallyto p 2 at small momentum. 
The width of the spectral line is instead controlled by the imaginary part of the self-energy, 
and at this order in 1/p, it scales like p 4 for small momentum: at this level of accuracy we 
recover Rayleigh scattering. The analytic calculation to this order can be confronted with 
numerical simulations (see section |T|), and it turns-out to be rather accurate for momenta 
of the order of the inverse interparticle distance. The density of states for small values of 
uj is also well reproduced. The maybe most interesting result of this paper is obtained in 
section [II C| , where the calculation of the self-energy is pushed up to order 1/p 2 . The order 



1/p 2 contribution to the spectral line- width is proportional to p 2 for small momentum, which 
is a qualitatively new feature. Even if the density of the system is high, for small-enough 
momenta the two loops contribution to the imaginary part of the self-energy is dominating. 
We believe this to be the underlying physical reason for the widespread experimental finding 
of p 2 scaling of the width of the dynamical structure factor. On the other hand, for larger 
momenta the order 1/p contribution is dominating, but it grows with p significantly more 
slowly than p 4 . This maybe of some importance in the experimentally relevant momentum 
range. 



A. The leading order 

Our strategy for computing the resolvent @ is straightforward: we consider it as the sum 
of a geometric series. One expands in series of 1/z the Fourier transform of the resolvent, 
see Eq.(0), in the following way: 

G(p,*) = -E(— ) R MW) (17) 

with 

MR (P) = ^ E (SkoM E - - /0*o - arfe) 

ko,ki. ..kit \ z± 
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) - f(x k 

ZR 



Since we are interested in the singularities (branch-cuts or poles) on the complex z plane of 
the resolvent G(p,z), see Eqs. ( p^3] , p~5D , it will be crucial to sum the series in Eq.fll8l) at all 



orders in 1/z. From now on we shall assume that the function / has a Fourier transform 
(denoted by / in the following) which is a rather severe limitation. In section [TV] we will 
see that, provided that particle correlations are taken into account on the calculation, this 
limitation should not bother us. One can readily check from Eq.flTgp that the contribution 
of the terms where the two particle labels in the argument of any / function are the same, 
vanish. This can be taken care of automatically by setting /(0) = 0. Since we have changed 
the value of the function at a single point, this will have no consequences while averaging 
over the positions of different particles. 

The very first step is to split the parenthesis in Eq . fll8f) , into 2 R sums (we keep the 
Kronecker deltas as a reminder of which hi index should be equal to which): 

e *(*fco-»* H ) 

M (P) = 2^ [<Wi/Ofc - Xzo)] [SkttefiXk! -Xzi)} •• • 

ko,Zo,...,Z R -l 



■ ■ [^k R -i,k R f{ x k R -i x zr^ x , 



e ip(x kQ -x kR ) 

+ Yl AT l-f( X ko - X kl )}[6 klM f(x kl - X Z1 )} 

fcO,fcl,«l,...,Zii-l 



• • • [3k R -l,k R f{ X k R -l X ZU-l)\ + 
e ip(x ko -x kR ) 

+ J2 jf [-f(xk -x kl )][-f(x kl -x k2 )}...[-f(x kR _ 1 -x ZR _ 1 )] (19) 

k ,k 1 ,...,k R 

The disorder average in Eq. flT5|) is simply obtained by performing the integrals: 

MW) = J II ^ MR (P) (20) 

Let us estimate the order of magnitude of the generic term in Eq.(|T9D, after the disorder 
average is performed. If there are no index repetitions, we will have an overall V~ R factor 
from the disorder average, because although there are R + 1 independent particles the 
translational invariance allows one to eliminate one integral. We have N\/(N — R — 1)! ~ 
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N R + l ways of choosing the particles labelling, and an extra iV _1 factor from the definition of 
the resolvent (||). Therefore, each of the above 2 R sums will be proportional to p R (p = N/V 
is the particle density). When two particle label coincide, we will have a missing N factor 
from the sums, and another missing V~ l factor from the position average, and the sum will 
be of order p R ~ l . It is thus clear that the M R {p) will be a polynomial of order R in p, with 
no term of order p°. In fact, the lowest order in p arise when there are only two different 
particle labels on the sums of Eq. (|i~9|) , which is easily seen to be linear in p: 

m e ¥) = p r i { r\p) + p^t^Im + ... + p4 r \p). (21) 

Since M R (p = 0) = 0, for every particle configuration, all the coefficient x\ R \p) in Eq.(^I|) 
vanish at zero external momentum, p = 0. 

Our zero'^-order approximation will consist in keeping only the leading term 

MW)^p R l { R R \p), (22) 

the first 1/p correction will consist in setting 

WHp) « p R 7™ (p) + p R ~ l j£ (p) , (23) 

and so on. 

Let us now calculate X^(p). As for any multidimensional integral, choosing with care 
the order in which the integrations are performed can tremendously simplify the task. In 
the case of Eq.(^) the matrix product yields a fairly natural ordering: from left to right, 
and when we encounter a diagonal term (Sk it k i+1 ) we shall first average the position of the Z{ 
particle. One easily shows that (recall that the Kronecker delta is simply a reminder, since 
ki is actually equal to fcj+i): 

J dx Zi ^8 khki+ J{x ki - x Zi ) = e**«.i/(0) , (24) 
J dx ki ^[-f(x ki -x ki+1 )} = e**H.i[-/(p)] • 

Therefore every integral simply shifts the external moment, p, one step forward, and leaves 
behind a /(0) or a —f(p)- When one arrives to the R+l particle, there is no integral left 
to be done, and both exponentials cancel. Then we find that 
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p R zg ) (p)=\pf(0)-pf(p) 



(25) 



and therefore we obtain the result already anticipated in Eq. (^3|) 



G (p,z) 



1 



(26) 



z - e(p) ' 



We see that at leading order, a plane wave with momentum p is actually an eigenstate of 
the matrix M with eigenvalue E, and the disorder does not play any relevant role. In other 
words, inside a wavelength 2ir/p there is always an infinite number of particles, smoothing 
out the density fluctuations of the particles: the system reacts as an elastic medium. 

Let us finally obtain the density of states at this level of accuracy, using Eq (|l2|) and G : 



for the dynamical structure factor: from the density of states one would say that the dis- 
persion relation is Einstein's like, without any momentum dependence! The way out of this 
contradiction is of course that in the limit of infinite p both e(p) and pf(0) diverge. The 
delta function in eq. (B7|) is the leading term in p, while the states which contribute to the 
dynamical structure factor appear only in the subleading terms in the density of states. 



For the calculation of the Zr_i(p) term, we need to consider only one particle label rep- 
etition. In other words, if we call Nr the average number of particles within the interaction 
range of the potential, we are calculating the 1/Nr correction of the formula M R = M , 
which neglects the statistical fluctuation of the matrix M due to disorder. 

The easiest way to proceed is as follows: one first identifies the places where the particle 
label repetition arises in the sequence of particle labels 




(27) 



We obtain a single delta function at p/(0), which is somehow contradictory with our result 



B. One loop 



(...1. ..!...), 
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and then one shifts the moment using Eqs. fl24|) until the place of the first repeated particle 
label. Then, one applies the backward version of the moment shift idea, from until the 
second particle label repetition, and one finally integrates over the position of the repeated 
particle. In this way one gets: 



P^'iflM _ ^ M/(o) - KpW W(o) - KpW 

z c+l 

l2 K/(0)-/(ff))]' 



(28) 

41 a+b+c=R-2 ^ * 

\ I ! < I I 1 I I I 1/(1 r ^ 

-«)-/(«))]' zi+1 , 

which is easily seen to yield 

£ pR '^B M = t G °(^)] 2 - n I iih G ^) Wip - ?) - fm 2 ■ (29) 

R=Q Z p J {Z7T) 

We want to obtain the self-energy term in (|T5|), which will receive contributions of all 
powers in 1/p: 

E{ P ,z) = E 1 (p,z) + E 2 (p,z) + ... (30) 

The relation between the above expresion and the resolvent is given by the Dyson resum- 
mation: 

G(p, z) - G Q (p, z) = G (p, z) Ei(p, z) + G (p, z)El(p, z) + E 2 (p, «) + ... G (p, 2) • (31) 
Comparing Eqs.(^TJ) and ( p9f) , we conclude that 

d. 



MP,z) = -J^ d Go(q,z)[pf(p-q)-pf(q)] 2 . (32) 

The 1/p 2 contributions to the self-energy will provide the Go(p, z )E 2 (p, z) term of Eq.([3ll), 
and also the more physically interesting £2(2?, z) term. 

Let us study in details the low exchanged momentum limit of Eq.([32|). It is clear that 
at p = the self-energy vanishes, as required by the translational invariance. We need to 
expand f(p — q) for small p, which due to the spherical symmetry of / yields 

f(p ~q) = /(g) ~(p-q)^f- + 0(p 2 ) , (33) 

= f{q) + {p-q) — + 0{p 2 ) , (34) 
qp 
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where we have used Eg . (|14D for the derivative of /. Substituting (|34|) in (p2|), and performing 
explicitly the trivial angular integrations in dimensions d we obtain 

MP,z)^p 2 —^—- r^ lJ ^C> (35) 



pd7i d / 2 T(d/2) Jo Z-e(q) 
, 2 1 - d /*(«=«>) [g(e)] d ~ 
' pd7r d /*T(d/2) Jo e q'{e){z-e)' 



t-J^f^*MtL.. (36) 

(37) 



In the last eguation, we have denoted with g(e) the inverse of the function e(q). Setting now 
z = E + i0 + , and observing that e(p) ~ Ap 2 for small p, we readily obtain 

ReE 1 (p,E + iO + )^p 2 —-=—rn^P ^^TTtI T> ( 38 ) 



pd7r d / 2 T(d/2) Jo q'(e)(E-e) 
m , E + ^ )x - -^l^- ) Aq( E)Y. (39) 

Since the principal part is a number of order one, the real part of the self-energy scales 
like p 2 (possibly with logarithmic corrections), and thus the speed of sound of the system 
renormalizes due to the 1/ p corrections. As a conseguence, the function q(E) is proportional 
to E 1 ! 2 ~ p at the maximum of the function of p Se(p, E), and the width of the peak of the 
Se(p,E) will scale like p d+2 . It is then easy to check (see (g)) that in freguency space the 
width of the spectral line will scale like 

T oc p d+1 , (40) 

as one would expect from Rayleigh scattering considerations. The result (f40"D for the asymp- 
totic regime p << 1 has been found at the one loop level. In order to predict correctly the 
spectral properties at very low external momentum p, it turns out that one must study the 
behaviour of the two loop contribution, as we shall see in the next section. Nevertheless, 
the one loop result is already a good starting point to perform detailed comparisons with 
the numerical simulations. To this end it will be useful in the following to introduce a 
concrete example with a specific choice of the function f(x), allowing a complete analytical 
computation. 
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1. The Gaussian case 



Our choice for the function / shall be the following 



■2 



fix) = e~*? , (41) 
f{p) = (2va 2 ) d l 2 e-^ . (42) 

The parameter a sets our length-scale, and for phenomenological purposes might be identi- 
fied with the inverse of the first maximum of the static structure factor. With this choice, 
we have 



i max = lim e(q) = p(2na 2 ) d/2 , 

q— >oo 

e(q)=e max fl-e"^) , (43) 



It is then straightforward to obtain in three dimensions 



Elfaz) = ^r- d£ £W4p), ( 44) 

\ n Jo z — e 



Q{q,p) = qe a 



22! 
2 



x + c _ a 2 p 2 sh(2a 2 pg) _ 2e _ CT 2Z sh(a 2 pq) 



2a 2 pq a 2 pq 
Therefore, setting z = E + i0 + , we find 



^fes+ioW^W'™"*^^ (45) 

V 7T Jo E — e 



/mE 1 (p,£ + iO + ) = —V27Tcr Q(q(E),p) . (46) 

One can easily check that the imaginary part of the self-energy on the peak is of order p d+2 , 
as previously announced. 

Turning now to the density of states, using Eq.([L^), Eqs.(^) and fllBj ), simplify now to 



/2a 2 r e ma X n(e) e ~ a q ^ ' 2 

ReE 1 (oo,E + iO + ) = J — P / rfe^| , (47) 

V 7r Jo E — e 

JmEi(oo, -E + i0+) = -v / 2Ta 2 g( J E) e^* 2 ^/ 2 . (48) 
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The zero" 1 order approximation, Qe^E) = S(E — e max ), Eq.(^7]), is largely modified at one 
loop. The imaginary part of the self-energy vanishes at and e max since both Q(q(0),p) and 
Q{<l{^max)iP) are zero. The zero th order delta function is moved to a value E* which verifies 
the relation^: 

E* =e max + ReZ 1 (^,E*), (49) 
and the density of states at this order is 



'2a 2 q(E) e" CT ^ I 2 

~ (E - ^(oo , E)f + q(E) 2 e~^ E ) 2 



9e(E) = - - f ^Z 2 ■ _ r ^__^ im2 , < E < e max , (50) 



9e{E) = ^ * i S(E -E*),E> e max . (51) 



£Z[i?e£i(oo,£)] 




dE 


E=E* 



C. Two loops 

The order 1/p 2 arises from including 2^ 2 ip) on the geometric series for the resolvent. 
The particle label repetitions can basically arise in four ways, that we schematically depict 

as 

(52) 
(53) 
(54) 
(55) 



1 


?, 


? 


1 


1 


?, 


1 


? 


1 


1 


1 




1 


1 


?, 


?, 



The repetitions depicted in ( [52] ) and ( |53| ) are genuine contributions to the self-energy. The 
repetitions ( |55| ) are contributions to the Dyson resumation of Si, (the Gq{j>, z)T,l(p, z) term 
in Eq.(pT|)), while the repetitions fl54]) contribute both to S 2 and to Gl(p, z)£f(p, z). 

In Appendix A we list all the two-loops terms giving the contribution to the self energy 
£2. Because of the complexity of the result, we shall focus herafter only onto the low 
momentum regime. For p = 0, there is a cancellation among all the various diagrams and 
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the self energy vanishes as it should. The leading contribution to the imaginary part of the 
self-energy scales as 

Jm£ 2 (p, E) = Cp 2 q d ~ 2 {E) (56) 

In the low energy, low momentum region, close to the peak p ~ E 1 ! 2 of the structure factor, 
it becomes: 

S 2 oc p 2 E^/ 2 oc p d . (57) 

Consequently, for small values of the external momentum p, the contribution to the width 
arising from the 1/p 2 order is much broader than the contribution from the 1/p order. In 
the frequency domain, the expression (|57D gives indeed the broadening: 

T oc p d - x (58) 

i.e. p 2 in the three dimensional case, very different from the result (fiO|). 

Putting all together, our final expression for the imaginary part of the self-energy in the 
low momentum regime, for a function fir) = f(r/a), is up to second order: 



JmS = pf(0) 



■ - d (pa) d+2 A, + (-U (pa ) d A 2 



(59) 



per" y pa u 

where the coefficients A\ and A 2 are pure numbers 0(1) which can be determined after 
choosing the adimensional function f(x) (the gaussian case will be studied in the next 
section). From the qualitative point of view it is worthwhile to note that (|59|) describes the 
most general case. It is quite easy to realize that even considering other multi-loops terms, 
for a small momentum, the leading order term would be of the form (EH). 

Let us note that the parameter of such expansions is actually the inverse of the number 
of particles = pa d , lying inside a cube whose dimensions are given by the range of the 
interaction. Let us picture the whole scenario as it emerges from two-loops computation: 

• As a very general feature of topologically disordered systems, we obtain the broadening 
oc p d ~ l of the resonance peak of the S(p, u) at low enough momentum p, where the two- 
loop contributions is dominant. Because of the behaviour, topologically disordered 
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systems spread the energy of an incoming plane wave on a wider range of frequencies 
than non-topologically disordered systems. This is a quantitative characterization of 
the intuitive impression that they are "more disordered". However, this feature is 
rather subtle, and only appears when the contributions oc to the self energy are 

considered. 

• At very high densities there exists a characteristic values p c where there is a cross-over 
from the p d ~ x regime to the oc p d+1 regime. From fl59|) the cross-over momentum is 
easily found: 

• For p > p c the one loop term mainly contributes and the broadening is controlled by 
Si. In this regime, the asymptotic result ( flO|) has to be taken with a lot of care. In fact 
it is supposed to be correct only for op « 1, a regime which basically is dominated by 
the two-loops contributions. If ap c is not small enough to be still in the same regime 
we are no longer legitimate to assume that (f|0]) give the correct results and we have 
to consider the full one-loop result fl32|). 

III. COMPARISON WITH NUMERICAL SIMULATIONS 

The dynamical structure factor Se{q, E) can be numerically computed. Its definition, 
Eq.(P), implies that it is a well normalized distribution function, 

/oo 
dES E (p,E) = l. (61) 
-oo 

The moments of this distribution function can be calculated in terms of the powers of the 
matrix M defined in Eq.([5|), 



/oo I 
dE S E (p, E)E R = -Y, eW**-"*> (M*) . j . (62) 
7V ij 
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If the function / is non-negative, one can cut the previous integrals at E = 0, since the 
quadratic form associated to the matrix M is semi positive-definite: 

Y^ViMijVi = ^Y,f( x i ~ - Vj) 2 > (63) 

Similarly, one can obtain the density of states replacing Eq.flB"!^) by 

I™ dE g E (E) E R = W ^(M«) y , (64) 

where the Vi are random numbers chosen with uniform probability between —1 and 1, and 
the overline now means average over the particle positions and the t>j. Thus we see that one 
can uses the method of moments^ complemented with a truncation procedure^. Using the 
gaussian function / of Eq.QHD truncated at four cr's, we have been able to reconstruct the 
Se(p, E) on systems with up to 32768 (32 3 ) particles, using 100 moments. All the simulations 
in this paper have been done generating ten samples of the disordered configurations: on 
a box of side L = 32a, we place at random pL 3 particles. We apply periodic boundary 
conditions, the minimum available momentum thus being J^. Regarding the choice of the 
density, let us recall that our natural length-unit is a, that can be considered as the analogous 
of the wavelength corresponding to the first maximum of the dynamical structure factor. It 
is therefore clear that the densities that correspond with the experimental situation are of 
the order of a~ 3 . In this paper we have explored the range 0.2a -3 < p < 1ct~ 3 . In what 
follows, momentum will be measured in units of a~ l and density in units of ct~ 3 . 

The choice of the number of calculated moments is conditioned by two conflicting goals. 
On the first place, the irregularities on the reconstructed probability distribution grows 
significantly with the number of calculated momenta, which is due either to the statistical 
fluctuations on the finite number of generated samples, or to round-off errors on the Lanczos 
recursion. On the other hand using a too small number of moments can smooth-out real 
features of the curves. As a rule we show the curve reconstructed with the minimum number 
of moments for which the width of the Se{p, E) is stable. We should emphasize that the 
moments method is a statistical one. With a finite number of particles, there are only pL 3 
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available eigenvalues. If the number of eigenvalues around the maximum of the Se{p,E) 
is small, the method of moments does not yield meaningful results. Moreover, since the 
Se{p,E) is an statistical quantity itself, meaningful results cannot be obtained even with a 
full diagonalization of the matrix M. In order to get a feeling on what difficulties might be 
encountered, it is useful to calculate the first and second moments of the S(p, E): 

dES E (p,E)E = e(p), (65) 
/_~ dE Se(p, E) E 2 = e(p) 2 + j^- d J dkf(k) if(k) - f(p - k)] . (66) 

It is clear in what sense in the p —>■ oo limit the Se(p, E) becomes a Dirac delta: its mean- 
value, e(p), grows like p, while its variance only grows like y/p. In particular, one can take 
the p — > oo limit in the above equations, and see what happens with the density of available 
eigenvalues. For the model described in Eq. (f|l]) , one obtains: 

\ dEg E (E)E = p(2na 2 f 2 , (67) 
Jo 

/•OO 

/ dEg E {E)E 2 = p 2 {2ito 2 Y + 2p{<na 2 fl 2 . (68) 
Jo 

The above results are rather discouraging, because one immediately observes that at fixed p, 
the larger is p, the lesser states are around e(p). Figure [I] may help in clarifying this point. 



19 



0.0 0.5 1.0 1.5 2.0 2.5 

- i i i i i I i i i i I i i i i I i i i i I i i i i i 




/ max 

FIG. 1. Density of states at densities p = 1 and 0.2. They have been obtained computing 
numerically 10 moments (see text). The X axis have been rescaled by e ma x, and the Y axis have 
been also rescaled in order to have unit area under the curves. 

The statistical sampling improves significantly with decreasing densities, where however 
our analytical calculation is less likely to work. In practice, the above considerations imply 
that the smaller is p, the smaller is the minimum momentum for which the Se{p, E) can be 
safely estimated on a finite box. 

As a first comparison between the analytical and the numerical computations, in figure 
(§) we show the DOS obtained numerically and the (continuous part) of the DOS at one- 
loop level as obtained by means of the eq. fl5T|) (neglecting the contribution of the delta 
function). As it has been already pointed out mill, the high density expansion is not really 
suitable to describe the DOS in the neighborhood of the maximum. In fact, we have shown 
that at one loop the zeroth order delta function splits from the continuous part of the 
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DOS, carrying still a finite weight. Furthermore the continuous part ( |51|) is defined only 
for E < e max ; in order to correctly compute the density of states in this region we need at 
least a partial resummation of the higher orders in the p l expansion, in the same way that 
it was done in for non translational invariant potentials. This task is not out of reach, 
however it goes beyond the scope of this paper. Because of this limitations, here only the 
low eigenvalues region of the spectrum can be safely compared. We see that the numerical 
curve obtained with a large number of moments is badly oscillating, because the single delta 
functions which form the spectrum of a finite matrix are correctly reproduced by the method 
of moments. If we consider the curve obtained considering only 10 moments instead, there 
is a smoothing of the previous curves and we see that it is in very nice agreement with the 
one loop computation. 

5 10 15 
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FIG. 2. Numerical DOS compared with the analytical prediction at one-loop level (neglecting 
the contribution of the delta function). The upper figure shows the comparison on the whole range 
of eigenvalues, and in the theoretical curve we have not drawn the delta function. The lower figure 
focuses onto the low eigenvalues part of the spectrum, showing the dependence of the numerical 
data with the number of moments. 

Let us turn to the comparison of the dynamical structure factor. One would like to 
check wether all the scaling regimes identified in subsection [11 (J| really hold. One should 



therefore know at which value of the exchanged momentum the one loop result will start to 
be dominating. The two coefficients of Eq.(^), for the Gaussian choice (|4l|) , in the 3 — d 
case the coefficients turn out to be: 

A x = 7/12fi 3 (69) 
A 2 & 1.25 0,1 (70) 

with 

The crossover momentum can then be estimated to be 0.31/cr for p = lcr~ 3 and 0.71/cr 
for p = 0.2cr~ 3 . In a box of size L = 32a the first available momentum is ~ 0.19/a, but 
unfortunately this is not the first momentum for which one can reconstruct accurately the 
Se(p,E). Therefore, we should expect that our data can be accurately described using 
only by the full one-loop expressions. Of course, one could try to compare with still lower 
densities, but this would be both risky (ours is a 1/p expansion) and unphysical (glasses have 
densities of order one, if measured in the natural length units of their interaction potentials). 
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FIG. 3. Dynamical structure factors at density p = 1, the momentum is given by p = 0.19n/cr, 
with n = 5,7, 9, 11. The numerical curves are obtained with 30 moments. 

The first comparison is for the density p — 1. As can be seen from Fig. [|, for a box with 
32 3 particles, when the position of the peak is at E peak < 1 — 2 there are too few states and 
the numerical results have no statistical significance. Hence, the lowest value of the external 
momentum which is possible to compare safely with the theory turns out to be p = 0.95/a 
(i.e. n = 5). The corresponding results are displayed in Fig. |^. We can see that for n = 5, 7 
the one-loop computation is in very good agreement with the numerical data. When the 
momentum is too large (n = 9,11 in the figure), E peak becomes comparable to e max and 
the agreement cannot be longer so good. The analytical curve indeed is defined only up to 
Cmax (at all orders in perturbation theory) and the tail which follows the peak cannot be 
reproduced by the theory without further developments. 

The next comparison is at density p = 0.2. Even if one might expect that this density 
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is too low for the 1/p expansion to be accurate, the one-loop computation still describes 
quite well the numerical curves, up to the region where E pea k and e max become comparable. 
One is tempted to conclude that the difference in fig. (|^) between theory and numerics for 
the curves which are not strongly affected by the tail effect (n = 4, 6) are due mainly to 
the following terms of the series in the 1/p expansion. In particular the correction of the 
broadening is of order ~ 10% while the position of the peak seems to be less affected by 
higher order corrections, specially for the n = 4 curve. 




1 2 3 4 5 6 



E 

FIG. 4. Dynamical structure factors at density p = 0.2, the momentum is given by p = 0.19n/<r, 
with n = 4, 6, 8, 10. The numerical curves are obtained with 10 moments. 

As we have pointed out in the introduction, a controversial point is the existence of a 
scaling law for the broadening of the spectral peak with respect to p, in the range of momenta 
explored by X-rays scattering experiments. On the other hand, there is a general agreement 
on the existence of the dispersion law for the peak position of the kind E pea k oc p 2 [uj peak cx p), 
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up to momenta where the corresponding wavelength is comparable with the mean separation 
of the particles, which enable us to identify the peak as an acoustic one. Having shown that 
our numerical data are reasonably described by the one loop computation, let us see what 
are the predictions for these spectral features at this level of accuracy. In the figure (|^) we 
plot the real and the imaginary part of Si at the point E pea k, obtained by means of (EHHBf) 
as functions of p, for the density p = 1 

0.05 0.10 0.50 1.00 5.00 



CD 



10° 




1 1 1 1 1 1 1 1 1 1 1 


iii i 
• 


i 












10 _1 






o 




10~ 2 






• 1 loop 




10~ 3 




p 

i i 1 i i i i 1 i i i 


O numerical 

iii i 


i 



10 l 



10 
10" 
10" 

10~ 



-2 



V 



• 1 loop 

O numerical 

■ 1 loop + 2 loop asymp. 



0.05 0.10 



0.50 



1.00 



5.00 



per 



FIG. 5. The real and the imaginary part of the self-energy computed at E pea k, for density p = 1. 
The open circles are the numerical results, which cannot be extended safely at lowest momenta. 
Full circles are the results of the exact one loop computation. The squares contain the exact one 
loop computation corrected by the small p expansion of the two loop result. 



The real part follows the dispersion law oc p 2 up to momentum of order 1/cr, and no 
significant dependence on the density is observed while that scaling law still holds, see 
Eqs.(|38|) and (p9|). For the imaginary part, we already know from the asymptotic analysis 
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of the p — > regime, that the one-loop low momentum behaviour is p 5 (p 4 in the frequency 
domain), however the range of momenta which are studied by numerical simulations, (and 
by IXS scattering experiments!) surely are not in that asymptotic regime. The figure shows 
clearly that the lowest momentum that we are able to numerically study is already in a 
region where that scaling law does not hold anymore. As a matter of fact, the embarrassing 
conclusion is reached that the scaling of the width of the dynamical structure factor at one 
loop (which coincides with the numerical simulation) could mimick a p 2 behaviour if the 
range of momenta which is studied is not wide enough. 



IV. ABOUT PARTICLE CORRELATIONS 



The formula (|TT| ) is a good starting point to speculate about the possibility of extending 
our computation to the case where the particles are not chosen with a random distribution 
but, for example, with the equilibrium Gibbs distribution. In that case, the average on the 
position of the particles is made usingS: 

i r R+l 

yR^v j lK'V// "-i.n..., •„.:). (72) 

i 

where the R + 1-points correlation function. Although the computation using the 

full correlation function would be exceedingly difficult, some progress can be made by using 
the so called superposition approximation: 

g(xt . . . x R+ i) = g(x 1 - x 2 )g(x 2 - x 3 ) . . . g(x R - x R+1 ) (73) 

where the pair correlation function is used to take into account the correlation of the position 
of the particles. The superposition approximation has been probed to be reasonably correct 
in describing the high order static correlation functions of the supercooled liquids, both 
in the computation of the coupling coefficients in the mode coupling theory^ and in the 
computation of the vibrational entropy@'@. 



26 



It is not difficult to see that, to the lowest order in the superposition approximation 
can be embebed in our calculation if we make the substitution: 



f(r) - g(r)f(r) (74) 

This is rather important, because, for typical applications, the function / being badly diver- 
gent at low distances, does not a have a Fourier transform. On the other hand, the function 
g(r) typically tends to zero at the origin exponentially, thus taking care of the algebraic 
divergence of /(r). 

However, the substitution (f74[) cannot be blindly applied to the calculation of the 1/p 
corrections. Let us consider for instance the second moment (|66|), that in the correlated case 
yield (let F be the Fourier transform of f(r)g(r)) 

e(p) = pF(0) - pF(p) , (75) 
Jp(pj = [e»] 2 + 2pJ dxg{x)f{x) (l - e vpx ) (76) 

while the subsitution (171) would have produced a factor g 2 (x)f 2 (x) in the above integral. 
A little thought reveals that, for the 1/p calculation reported in Eq. (Effi) , this problem will 
arise whenever there is a string of only diagonal terms ($ki,k i+1 f( x ki ~ x zi)) between the two 
repeated particle labels. The diagonal terms are simply [pF(0)] b / z b+1 that add up to 

1 

z - pF(0) ' 

Therefore, we conclude that the self-energy at order 1/p in the superposition approximation 



is 



s > te z) = \ I - " - " FM ' 2 (t^W - T^m) + (77) 



+ 2p / dxg{x)f (x 



, ; 1 1 - e ipx ) 



z - pF(0) ' 

It is clear that none of the two new terms contribute to the imaginary part close to e(p), 
and therefore our discussion of the broadening of the spectral line, at this order of the 1/p 
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expansion is basically unchanged. One could similarly discuss the peculiarities of the 1/p 2 
term, but this is left for future work. However, the single term 

2pj dxg(x)f(x) [ z _ pF{ > y (78) 

points to the possibility of a qualitatively new behaviour. Indeed, this term is the first of a 
full category, in which the repeated particle labels always corresponds to the same pair of 
particles, separated by a string of diagonal terms. It is not hard to show that the contribution 
to the self-energy of this family of terms is: 

fdx g(x) ( 1 - e ipx ) } 2f ^XZ =2p fdx g(x)f 2 (x) - . . (79) 

This result implies that when the function f(x) diverge at small distance, the density of state 
will have a tail to infinity, exponentially supressed by the g(x) function. This exponential 
tail is of a different origin from the instantonic contribution recently calculated in Ref.0. 
Indeed the instantonic exponential tail arise also for a fully regular f(x), due to anomalously 
dense regions. 

V. DISCUSSION, CONCLUSIONS AND OUTLOOKS 

The aim of this work was to study the spectral properties of euclidean random matrices 
using an expansion on the inverse particle density number. That approach is supposed 
to describe correctly topologically disordered systems, in the context of the Instantaneous 
Normal Modes approach to their high-frequency dynamics. In particular, we have performed 
a two-loop computation of the dynamical structure factor, where the expansion parameter 
is the inverse of the density of particles. The results have been analytically studied in the 
low momentum regime, where the disorder reveals itself by the broadening of the resonance 
peak. The width of that spectral line has been computed up to order 1/iVjj, where Nr is the 
number of particles which effectively interact with a given particle. The behaviour of the 
width r in the frequency domain, as function of the momentum p has been shown to follow a 
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quite complicated law. When the exchanged momentum, p, is quite smaller that the inverse 
wavelength of the first maximum of the static structure factor, the width of the spectral line 
follows r oc p 2 . This seem to be an intrinsic difference between topologically-ordered and 
topologically-disordered random systems, the Rayleigh p A scaling being always found for the 
formers. At a momentum depending on the particle density and the details of the potential, 
the 1/Nr contribution to the self-energy becomes dominant, and the p A scaling appears. 
Finally, when the exchanged momentum is of the order of the inverse wavelength of the first 
maximum of the static structure factor, the dynamical structure factor starts its collapse 
onto the density of states. The width of the spectral line cannot grow at such a violent 
pace as p A , and if a not too large momentum-range is explored, it can certainly mimick a p 2 
scaling. It is possible that the p 2 scaling found in the IXS experiments correspond to the later 
transient behaviour. Nevertheless, it would be very interesting to explore experimentally a 
wider range of exchanged momentum (maybe combining several scattering probes), in order 
to test this scaling picture with three separated regimes. 

It would be also very useful to compare the predictions of this computation with the 
ones which have been obtained in the framework of mode coupling approximation^. There 
an hard sphere system is investigated and a p 2 law is obtained only for very low values of 
the frequency of the peak, which at the moment are not accessible to the experiments, while 
a crossover to a quite different behaviour is present when the external momentum moves to 
higher values, comparable with the ones explored by the IXS scattering experiments. Even 
in that computation no clear scaling law actually emerges and qualitatively that scenario is 
quite similar to the one emerging from our approach. A quantitative comparison though, 
is not straightforward because our analytic computation does not concern any particular 



model of topologically disordered system, the choice ([11]) being quite generic in order to 
understand which are the main features which determine the particular kind of broadening 
of the resonance peak. 

The main point, however, is the reliability of our computation in understanding the 
experimental findings and, possibly, in predicting features not yet observed. Apparently, if 
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one wants to apply the computation presented here to the study of the spectral properties 
of realistic systems, such as silica or fragile glasses, one faces the following limitations: 

• The computation has been explicitly performed in the case where particles are uncor- 
rected which is deeply different from the situation found in glasses and supercooled 
liquids. 

• The values of the coefficients A± and A<i in ([59]) seems to depend very badly on the 
choice of the interaction f(r) hence at this level is quite unclear the generality of 
the predictions that one could obtain, from (|60|) before performing the appropriate 
computations for a given system. 



However, in section [TV] it has been shown that the correlations between the particles can be 
at least partially taken into account by performing the transformation (|74]), which amounts 
to dress the bare interaction f(r) between two particles at a given distance r, with the 
probability, given by the static structure factor g{r) of the real system to be studied, to 
find them at that distance r. Furthermore this procedure corresponds to the well-known 
superposition approximation of the correlation functions of the theory of liquids, a kind of 
approximation that has been largely exploited previously in different analytical approaches 
to the study of the glassy phase. More interestingly, the superposition approximation makes 
even the latter limitation less severe. As a matter of fact the Gaussian choice ([|]]) has 
to mimic not the bare interaction f(r), that would be an hopeless task, but the function 
f(r)g(r) which sounds more feasible. We have also shown how the terms for which the 
substitution /(r) — > f(r)g(r) cannot be made, give rise to an exponential tail of the density 
of states, of non instantonic (seeil) origin. More in detail, by considering only the bare 



interaction one should take: 



where V(r) is the pair potential of the particles, which typically has a strong repulsive 
barrier at r — * 0, missing in the Gaussian choice. The static structure factor g(r), on the 
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other hand, is ~ e ^( r ) at r — > 0, hence the function f(r)g(r) is typically formed by a single 
peak in the position r , where the static structure factor has its maximum, is not singular 



at r — * and it does have a finite Fourier transform looking quite similar to (0). As a 
first approximation, the above argument yields a ~ ro, i.e. a 'few' Angstrom for realistic 
systems, implying that X-rays experiments surely cannot explore the low momentum region 
but they can observe the region pa > 0.2 — 0.3, which is roughly the same region as that 
investigated in the numerical simulations. 

Hence we expect that the computation we have presented here, if suitably modified 
for including the details of the system under study, will allow to compute the dynamical 
structure factor of structural glasses and supercooled liquids. 
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APPENDIX A 



In this appendix we shall show the two loops integrals and the corresponding diagrams 
arising from the 1/p expansion of the resolvent to order 1/p 2 . It turns out to be useful to 
express the integrals in terms of the propagators pf(q) and G s (q), where 
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Diagrams: C 
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Diagrams: E 
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Diagrams: H 

k 




^j d ^0GM)G s (k)G s { P - q -k) 

i / ji^ Gsiq) pf ~ {k) pf ~( p ~ q ~ k) 
^0gm) (pf(k)) 2 (pf( P - q -k)r 



Diagrams: J 
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Diagrams: K 
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